# SETWD ===========================
rm(list=ls())


# LIBRARIES ===========================

library(tidyverse)
library(foreign)
library(stargazer)
library(lmtest)
library(multiwayvcov)

# IMPORT DATA ===========================

dfgi <- read.csv("data/sve-08_04_24.csv",stringsAsFactors = F)

## TAB A9. Logit: Region Fixed Effects ----

temp <- dfgi%>%
  filter(!is.na(gini_disp))

m1 <- glm(erosion.strict ~ gini_disp + region, data=temp, family=binomial)
crse1 <- coeftest(m1, cluster.vcov(m1, temp$country.name))
crse1

m2 <- glm(erosion.strict ~ gini_disp + log(gdppc) + region, data=temp, family=binomial)
crse2 <- coeftest(m2, cluster.vcov(m2, temp$country.name))
crse2

m3 <- glm(erosion.strict ~ gini_disp + log(gdppc) + year + region, 
          data=temp, family=binomial)
crse3 <- coeftest(m3, cluster.vcov(m3, temp$country.name))
crse3


stargazer(crse1,crse2,crse3, title="Logit: Region FE",digits=3,
          header=F,df=F,omit.stat=c("f","ser","aic","ll"),
          star.cutoffs=c(0.05,0.01,0.001),
          add.lines = list(c("Observations", nobs(crse1), nobs(crse2), nobs(crse3))),
          star.char=c("*","**","***"),
          notes="$^{*} p<0.05$; $^{**} p<0.01$; $^{***} p<0.001$",
          covariate.labels=c("Gini","Logged GDP per capita",
                             "Year",
                             "Europe/Central Asia",
                             "LatAm/Caribbean",
                             "MENA",
                             "North America",
                             "South Asia",
                             "Sub-Saharan Africa"),
          notes.append=F,table.placement="h!",
          out="tables/apptab-regionFE-07_04_23.tex")


## TAB A10. Logit: All Controls ----

temp <- dfgi%>%
  filter(!is.na(gini_disp),!is.na(gdppc),!is.na(prs.ge),!is.na(v2cacamps),
         !is.na(democracy_duration))

m6 <- glm(erosion.strict ~ gini_disp + log(gdppc) + year + prs.ge + 
            democracy_duration + v2cacamps + region, 
          data=temp, family=binomial)
crse6 <- coeftest(m6, cluster.vcov(m6, temp$country.name))
crse6

stargazer(crse6, title="Logit: All Controls",digits=3,
          header=F,df=F,omit.stat=c("f","ser","aic","ll"),
          add.lines = list(c("Observations", nobs(crse6))),
          star.cutoffs=c(0.05,0.01,0.001),
          star.char=c("*","**","***"),
          notes="$^{*} p<0.05$; $^{**} p<0.01$; $^{***} p<0.001$",
          covariate.labels=c("Gini","Logged GDP per capita",
                             "Year",
                             "Bureaucratic Quality",
                             "Age of Democracy",
                             "Political Polarization",
                             "Europe/Central Asia",
                             "LatAm/Caribbean",
                             "MENA",
                             "North America",
                             "South Asia",
                             "Sub-Saharan Africa"),
          notes.append=F,table.placement="h!",
          out="tables/apptab-allcontrols-10_27_24.tex")


## TAB A11. Logit: Alternative Demonstration Estimates ----

erosion.cdf <- dfgi%>%
  group_by(year)%>%
  summarise(erosion=sum(erosion.strict,na.rm=T))%>%
  mutate(erosioncdf=cumsum(erosion))%>%
  dplyr::select(year,erosioncdf)%>%
  mutate(year=year+1)

dfgi2 <- dfgi%>%
  left_join(erosion.cdf)

dfgi2$erosioncdf[dfgi2$year==1995] <- 0

temp <- dfgi2%>%
  filter(!is.na(gini_disp),!is.na(gdppc))

m1 <- glm(erosion.strict ~ gini_disp + log(gdppc) + erosioncdf, 
          data=temp, family=binomial)
crse1 <- coeftest(m1, cluster.vcov(m1, temp$country.name))
crse1
nobs(crse1)



m2 <- glm(erosion.strict ~ gini_disp + log(gdppc) + as.factor(year), 
          data=temp, family=binomial)
crse2 <- coeftest(m2, cluster.vcov(m2, temp$country.name))
crse2
nobs(crse2)


stargazer(crse1, crse2, title="Logit: Alternative Demonstration Estimates",digits=3,
          header=F,df=F,omit.stat=c("f","ser","aic","ll"),
          add.lines = list(c("Observations", nobs(crse1), nobs(crse2))),
          star.cutoffs=c(0.05,0.01,0.001),
          star.char=c("*","**","***"),
          notes="$^{*} p<0.05$; $^{**} p<0.01$; $^{***} p<0.001$",
          covariate.labels=c("Gini","Logged GDP per capita",
                             "Cum. Prior Erosion Years"),
          notes.append=F,table.placement="h!",
          out="tables/apptab-yearFE-08_05_24.tex")



